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ABSTRACT 

Wave propagation in inhomogeneous media has been studied for such diverse applications as propagation of 
radiowaves in atmosphere, light propagation through thin films and in inhomogeneous waveguides, flow visualization, and 
others. In recent years an increased interest has been developed in wave propagation through shocks in supersonic flows. 
Results of experiments conducted in the past few years has shown such interesting phenomena as a laser beam splitting and 
spreading. 

The paper describes a model constructed to propagate a laser beam through shock-like inhomogeneous media. 
Numerical techniques are presented to compute the beam through such media. The results of computation are presented, 
discussed, and compared with experimental data. 


1. INTRODUCTION 

Interest in a medium with a rapidly varying refractive index has been increasing recently partially due to the advent 
of supersonic flight, a growing need for better flow visualization systems and a deeper understanding of light propagation 
through shocks. In that respect, attempts have been made to explain the refraction phenomenon 1 and formation of refractive 
fringes 2 and to conduct mathematical and experimental analysis of light diffraction on and transmission through plane shock 
waves 3,4 . Such phenomena as light diffraction 5 and scattering 6 ’ 7 on shocks have been observed and reported. Also, 
experiments have been performed to determine a normal shock location. 8 9 In view of this development in the experimental 
field of shocks visualization and analysis a need has arisen for a deeper understanding of the phenomenon of light propagation 
through a highly inhomogeneous medium. Thus, theoretical and computational models to perform numerical analysis have 
become important for explaining the recently observed phenomena such as laser beam splitting and broadening. 

The purpose of this paper is to present a computational model of a laser beam striking an inhomogeneous body under 
a grazing incidence. The model includes the inhomogeneous body, the incident laser beam, and a computational scheme to 
propagate the beam through the inhomogeneity under the grazing incidence. 
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2. DESCRIPTION OF THE MODEL 


To evaluate the phenomenon of wave propagation through inhomogeneous media the following model has been 
constructed. The inhomogeneous media is assumed to be a penetrable circular cylinder with a cylindrically symmetric profile 
of the refractive index. The radial distribution of the refractive index profile has a shock- like profile. Such a profile has been 
described in the literature: 10 13 


where 


n < r > = n low + ■ 


An 


1 + exp\ - 4 


r-R 


** = n high - n io» . 

n high n iow are the maximum and minimum values of the refractive index respectively, 

R is the radius of the inhomogeneous cylinder, 

__ — 

r = y x + y is the radial coordinate, 

jc and y are Cartesian coordinates of the point of observation, 

L is the shock thickness. 


( 1 ) 


Parameters An , R , and L describe the shock-like profile of the refractive index. Figure 1 represents an 
example of a 2-dimensional distribution of the refractive index with An = 0.01 , R = 25A 0 , and L = A 0 , where A 0 is the 

wavelength in vacuum. In this work we assume that n hw = 1 and L = 0. Thus, when L = 0, we have a homogeneous 

cylinder with the index of refraction n = l + An placed in another homogeneous medium with the index of refraction equal 
to 1. 


The cylinder described above is placed in the Cartesian coordinate system with its long axis along the vertical Z axis. It is 
illuminated by an incident electromagnetic field with the propagation vector normal to the long axis of the cylinder. The 
electromagnetic field is a sheet of light with a constant intensity in the direction along the axis of the cylinder and the 
Gaussian intensity profile in the direction normal to it. We will call this sheet of light a laser beam. The electic and magnetic 
field vectors are chosen to form a transverse magnetic (TM) wave. Assuming that the direction of propagation is the Y 
direction, the intensity of the two dimensional incident field can be written as: 14,15 


E(x,y,t)=A. 
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u> 0 is the Gaussian beam waist, and A is the wavelength of radiation in the medium with refractive index n. Thus 
A = A 0 /n , where A 0 is the wavelength of radiation in vacuum. Coefficient A is a normalization constant. 


The beam and the cylinder are positioned in such a way that the optical axis of the beam strikes the cylinder at 
the grazing incidence at r = R . Selection of the described above configuration permits reduction of a three dimensional 
problem of wave propagation to a two dimensional one. 

Selection of a small diameter incident Gaussian beam striking a relatively large diameter scatterer gives an 
opportunity to separate a scattered field from the incident one. The idea has been implemented for curvature radii 
measurements. 16 In that experiment a laser beam impinging at a grazing incidence on a surface produced a diffraction edge 
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wave. A reduction of the laser beam diameter led to a separation of the edge wave from the incident beam. In the region 
where the two fields overlapped, diffraction fringes were observed. 

To compute propagation of the incident beam through the medium described a hybrid method has been selected. 
The method consists of two parts, propagation through the inhomogeneity and projection of the emerged field into the far 
field. The first part of the problem, propagation of an electromagnetic field through an inhomogeneity, is computed using 
the finite-difference time domain (FD-TD) method. The wavefront that emerges as result of calculations then is propagated 
to a remotely located screen using the Fresnel diffraction equation. 

3. FINITE-DIFFERENCE TIME-DOMAIN METHOD 


3.1 Introduction 

The FD-TD method of computing the electromagnetic wave propagation is based on a simultaneous solution of a 
system of the first order partial differential equations derived from Maxwell’s time dependent curl equations. 17 19 
Furthermore, the electric and magnetic field components are positioned in a specific manner described by the Yee algorithm. 
The algorithm permits solving for both electric and magnetic fields in dme and space rather than solving for the electric field 
along with a wave equation. Those electric and magnetic components are positioned in space in a specific interleaved way 
which permits a natural satisfaction of tangential field continuity conditions at the interfaces. Due to the fact that the process 
of solving partial differential equations in an unbounded domain using discrete techniques involves a truncation of the 
solution domain, an approximated boundary is introduced at a finite distance from a scatterer. This approximate boundary 
condition is also called an absorption boundary condition (ABC). The ABCs developed by Mur 21 are specially designed to be 
used with the FD-TD method. Simultaneous discretization in space and time domains requires temporal stability. The time 
domain discretization scheme is stable if the ratio of spatial segmentation distance to the time step size satisfies the Courant 
criterion. 22 A straight forward application of a cubical Yee cell in Cartesian coordinates to curved surfaces leads to stair- 
case approximation. The resultant stepped edge profile of the approximated surface generates an error. 23 One of the ways to 
minimize the problem is to make the cells small. 


3.2 Numerical implementation 

To compute wave propagation using the FD-TD method the scattered field formulation has been chosen. It is based 
on splitting the total field on the known incident and unknown scattered fields, performing the FD-TD computation of the 
scattered field, and adding the incident field to it to obtain the total field. For a 2D problem pertaining to the model described 
in the previous section such formulation in a case of TM wave leads to the following equations: 
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The FD-TD discretization process applied to equations derived from scattered field formulation of TM wave 
propagation gives the following: 

E" + \iJ) = E”(i,j) 


'H^\i + l/2,j)-H; + y 2 (i-X/2J) 


-H n x +vl (i, j + 1/2) + H n x +l/Z (/, j - 1/2) 


n+V 2, 


, _ Incident , . 

l-e r dE z (t,j) , 

+ £ A/ , 


£ r dr 

H x +V2 (i,j + l/2) = H?T V 2 (i,j+ i/2)-C H [E^(i,j + l)-E^(i,j)] , 
Hy +V2 (i + 1/2 , 7) = H^~ yl (i + 1/2 , j) + C„ [e” 0 + 1, J) - £" 0, y)] , 


( 8 ) 

(9) 


where coefficients associated with electric and magnetic fields are correspondingly C E = ^ and C H - * 

and A = Ax = Ay . For simplicity, notations for scattered fields in the equations above are omitted. The absorption boundary 
conditions used are the Mur’s conditions of the 2 nd order for the edges of the computational domain and of the 1 st order for the 
comers. 


3.3 Computational results 


The computational domain is selected to be 60 wavelengths wide in X direction and 80 wavelengths long in the 
direction of beam propagation, or Y direction. To minimize a negative effect of the stair case approximation the size of space 
steps is chosen to be 0. 1 of the wavelength. This resulted in a two dimensional grid with 600 x 800 grid points. 2000 time 
steps are used to achieve a steady state of the computed field. The time step is selected according to the Courant criterion to 
be equal to 0.99 of the Courant number A/ c . The Courant number for a two-dimensional problem is derived from the Courant 
stability criterion : 


At 


c 



( 10 ) 


where 

A = Ax = Ay is the distance between the grid points, 
c Q = 1 j ^ je jx is the speed light in vacuum. 


£ 0 and /x 0 are permittivity and permeability of vacuum respectively. 


A two dimensional Gaussian beam having the wavelength A 0 = 1 fi and the waist radius w 0 = 10A 0 enters the 
computational domain located at a distance y 0 = 200A 0 from the waist. The cylinder and the Gaussian beam are positioned 
in the computational domain as presented in Figure 2. The propagation of the beam is in the Y direction in a such way that its 
directional axis passes through the middle of the computational domain. A cylindrical shape body with the refractive index 
different from the one of the surrounding medium is placed in the passage of the beam the way described in the previous 
sections. 

Figure 3 shows propagation of the Gaussian beam through a medium which contains a cylinder with the radius 
R = 30A 0 and the maximum refractive index difference between the cylinder and the surrounding medium An = 0.005 . 
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Splitting of the incident Gaussian beam and formation of a double peak and fringes are clearly seen on the picture. The 
phenomena are caused by a combination of effects. The most significant are the interference between the incident Gaussian 
beam and the diffracted edge wave and scattering by a dielectric cylinder. 

Computed effects of the radius of the cylinder and its refractive index on the wave propagation are presented in 
Figures 4 through 7. The first two of them, Figures 4 and 5, represent three-dimensional views of Gaussian beams propagated 
through cylinders similar to the one used to obtain data shown in Figure 3 but for An = 0.002 and An = 0.008 respectively. 
Calculated intensity distributions at the exit from the computational domain for An = 0.002, 0.005, and 0.008 are presented 
in Figure 6. Figure 7 shows the effect of a change in the radius R of the cylinder with An = 0.005 . 

Thus, the calculations have shown that the cylinder radius R and refractive index difference An have a significant 
effect on the relative amplitude of the two main peaks in the intensity distribution. It can be seen from the figures that an 
increase in any of these parameters leads to an increase in amplitudes of both peaks. Moreover, in response to changes in 
these parameters, the amplitude of the peak to the right changes more rapidly than the one to the left. Another factor that 
plays an important role in the intensity distribution is the relative position of the beam and the cylinder. 

4. FORMATION OF IMAGE IN THE FAR FIELD 


One of the methods to propagate optical fields involves the Fresnel diffraction integral. The integral facilitates 
propagation of an optical disturbance from one plane with coordinates £ and 7 ] to another one with coordinates x and y and 
located at a distance z from the first. Applying a conventional technique described by Weaver 24 to a two dimensional 
problem and maintaining the same coordinate notation as in the previous chapters, the following form of the Fresnel 
diffraction equation can be derived: 


V 2 = 


Ke yky 

V7 


where K is the inclination factor. 


Y 


(S)exp 




( 11 ) 


The last expression can be written in terms of the Fourier transform and then solved numerically. Using the 
established procedure the following is obtained: 

y 2 («) = e jky ■ 4>j («) • e ~ inXyul = 'Fj («) ■ H(u ) , 

where 


4^ (w) and ^(w) are the Fourier transforms of Y\( x ) an< ^ respectively, 


H(u) is the free space transfer function of the system, 


H(u) = e J J e 


jky -jnkyu 2 


( 12 ) 


Thus, the process of propagating an optical field from one location to another consists of three steps. These steps are 
computing the Fourier transform of the field in the original plane, multiplying it by the free space transfer function, and 
performing the inverse Fourier transformation of the resultant expression in order to find the field at a new location. To 
perform the direct and inverse Fourier transformations a Fast Fourier Transform algorithm, based on the Darnel son-Lanczos 
Lemma, and computer codes are adopted from available literature on numerical techniques. 25 Results of propagation are 
presented in the following figures. In the first series of figures the original field is computed using the FD-TD method and 
then propagated to distances of 20A 0 and 40A 0 . Figure 8 shows the intensity distribution at a distance of 20A 0 for cases 
when the refractive index difference An - 0.005 and An - 0.008 . The radius of the cylinder R in both cases remain the 
same, R = 30A 0 . When the distance increases the pattern goes through transformations. The sharp changes in the computed 
intensity distribution become smoother and eventually disappear. 

A phenomenon of beam spreading can be observed by comparing the intensity distributions of an undisturbed or 
reference Gaussian beam with the one that emerges after propagating through the cylinder. The beam spreading manifests in 
an increase of the spatial width of the curve that forms the intensity distribution. In Figure 9 the beam spreading can be seen 
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at the right side of the curve next to the reference Gaussian profile. The cylinder used to compute the data has the following 
parameters: R = 30A 0 and An = 0.008 . This phenomenon has already been observed experimentally and published. 9 Figures 

10 to 1 1 show intensity distributions at a distance of 8OA 0 from the exit from the computational domain from FD-TD 
computations with An = 0.005 and An = 0.008 respectively. The beam splitting and broadening are present. Fringes can also 
be seen. Increase in the refractive index and/or the radius of the cylinder will lead to enhancement of these phenomena. 

5. CONCLUSION 

A two-dimensional model and hybrid computational technique have been proposed in this paper to propagate a 
Gaussian beam through inhomogeneities with shock-like profiles into the far field under a grazing incident condition. 
Computing of the beam propagation through the computational domain is performed by the FD-TD method . The shape of 
inhomogeneity is selected to be cylindrical. The resultant fields are then propagated into the far field using the Fresnel 
diffraction equation and Fourier transformation. The computed patterns show effects of the refractive index and the radius 
of the cylinder. The patterns of intensity distribution of a Gaussian beam in the far field show beam splitting and spreading. 
These phenomena have been also observed experimentally. An example of the experimentally obtained intensity profile of a 
Gaussian beam after passing a bow shock is shown in Figure 12. The experimental setup used was similar to the one 
described in literature. 6,7 To generate bow shock a cylindrical blunt body was inserted in the supersonic flow. A laser beam 
sent through the shock under the grazing angle of incidence was projected to a remotely located screen. A CCD camera 
captured the image of the beam on the screen and displayed the beam intensity profile on a computer screen. The beam 
intensity profile clearly shows beam splitting and formation of fringes. Thus, the model and computational method are 
supported by experimental data. Moreover, the phenomenon of beam spreading by a shock may be used as the basis for 
shock detection. 

An extension of the method proposed in this paper into the three-dimensional domain will be one of the first future 
areas of effort. Inhomogeneous bodies then will be spheres with shock-like profiles of the refractive index and large 
diameters. To build a three dimensional computational model with geometrical dimensions close to those that appear under 
real conditions some shortcomings of the presented method have to be overcome. One of the shortcomings comes from the 
limitations of the FD-TD method. Methods based a phase object approximations 26,27 may help to eliminate those limitations. 
One of the methods, anomalous diffraction approximation, 28,29 is especially attractive when a phase object has it refractive 
index close to the one of the surrounding medium. However, the method would have to be modified to include potential 
refractive effects of the spheres. 

Other areas that will deserve future attention involve the large angle scattering and polarization phenomena. In 
order to increase the field of view and evaluate effects associated with large angle scatter, the inclination factor K has to be 
closely evaluated. The beam propagation into the far field using the Fresnel diffraction equation is based on a scalar field 
formulation. This means that polarization of the incident beam is not taken into account. Development of a vector field 
formulation and an associated computational technique represents a certain interest and challenge. 
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Figure 1 : Example of 2D distribution of the refractive index 
with An = 0.01 , R = 25A 0 , and L = A 0 . 


Figure 3: Results of computation of a Gaussian beam 

propagation through inhomogeneous media with 
R = 30A 0 and An = 0.005 ; grazing incidence 



Figure 2: Top view of the computational domain with 
relative orientation of the incident beam and 
cylinder; grazing incidence. 



Figure 4: Results of computation of a Gaussian beam 

propagation through inhomogeneous media with 
R = 30A 0 and An = 0.002 ; grazing incidence. 
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Relative intensity 



Figure 5: Results of computation of a Gaussian beam 

propagation through inhomogeneous media with 
R = 30A 0 ,and An = 0.008 ; grazing incidence. 



Distance in wavelengths, x 0.6 

Figure 7: Calculated intensity distributions at the exit from 
the computational domain for An = 0.005 and 
different radii R. 
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Figure 6: Calculated intensity distribution at the exit from 
the computational domain for An = 0.002 , 
0.005, and 0.008; ( R = 30A 0 ). 


Figure 8: Intensity distribution at 20A 0 distance for cases 
of the refractive index differences An = 0.005 
and An = 0.008 ; ( /? = 30A 0 ). 
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Distance in wavelengths, x 0.1 



Distance in wavelengths, x 0.1 


Figure 9: Beam spreading of a Gaussian beam at 20A 0 
distance by a dielectric cylinder ( R = 30A 0 
and An = 0.008 ) under grazing incidence. 


Figure 11: Intensity distribution at 80A 0 distance 

obtained using the FD-TD data and Fresnel 
diffraction equation; ( R = 30A 0 , An = 0.008 ). 



Distance in wavelengths, x 0 1 


Figure 10: Intensity distribution at 80A 0 distance 

obtained using the FD-TD data and Fresnel 
diffraction equation; ( R = 30A 0 , An = 0.005 ). 



Figure 12: Example of an experimentally obtained 
intensity profile of a Gaussian beam after 
passing through a bow shock. 
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